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Abstract 

In this paper we show how to compute algorithmically the full set of algebraically independent 
constraints for singular mechanical and field-theoretical models with polynomial Lagrangians. 
If a model under consideration is not singular as a whole but has domains of dynamical (field) 
variables where its Lagrangian becomes singular, then our approach allows to detect such do¬ 
mains and compute the relevant constraints. In doing so, we assume that the Lagrangian of a 
model is a differential polynomial and apply the differential Thomas decomposition algorithm to 
the Euler-Lagrange equations. 

Keywords: Singular Lagrangians, Lagrangian constraints, Euler-Lagrange equations, Thomas 
decomposition, differential algebra 


1. Introduction 

Models with singular Lagrangian^lay a fundamental role in modern elementary particle 
physics and quantum field theory (cf. IiIh]). Singularity of such models is caused by local sym¬ 
metries of their Lagrangians. Gauge symmetry is the most important type of local symmetries 
and it is imperative for all physical theories of fundamental interactions. The local symmetry 
transformations of differential equations describing dynamical systems or field-theoretical mod¬ 
els relate their solutions satisfying the same initial (Cauchy) data. Eor dynamical systems with 
only one independent variable the initial data include (generalized) coordinates and velocities 
whereas for field-theoretical models they include the field variables, their spatial and the first- 
order temporal derivatives (‘velocities’). The presence of local symmetries in a singular model 
implies that its general solution satisfying the initial data depends on arbitrary functions. 

A distinctive feature of singular Lagrangian models is that their dynamics is governed by 
the Euler-Lagrange equations which have differential consequences in the form of internal (hid¬ 
den) constraints for the initial data. This is in contrast to regular constrained dynamics whose 
constraints are external with respect to the Euler-Lagrange equations. 

Given a model Lagrangian, it is very important to verify whether it is singular, and if so 
to compute internal constraints that are hidden in the Euler-Lagrange equations. Knowledge of 
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such constraints is necessary in constrained dynamics for the local symmetry analysis |[ll] and for 
quantization ^ of the singular model. As it was observed first in Q, 0] in the framework of the 
Dirac formalism for investigation and quantization of constrained dynamical systems and gauge 
field theories, the Dirac procedure of the constraint calculation is nothing else than completion 
of the Hamiltonian equation system to involution. In this respect, the Hamiltonian constraints are 
integrability conditions |0. The algorithmic aspects of computation and separation of Hamilto¬ 
nian constraints into the first and second classes were studied in 10. Computation of independent 
Lagrangian constraints as integrability conditions for the Euler-Lagrange equations is essential 
for well-posedness of initial value problems lll0i . Such constraints are detected in the course 
of completion of the Euler-Lagrange equations to involution. 

In the present paper we consider Lagrangian models whose Lagrangians (mechanics) and 
Lagrangian densities (field theo^) are differential polynomials. Under this condition we show 
that the Thomas decomposition |j3[I3]j being a characteristic one (cf. d) for the radical differ¬ 
ential ideal generated by the polynomials in Euler-Lagrange equations, provides an algorithmic 
tool for the computation of Lagrangian constraints. Unlike the traditional linear algebra based 
methodology used in theoretical and mathematical physics fll] for computing linearly indepen¬ 
dent Lagrangian constraints, our approach is not only fully algorithmic but also takes into account 
rank dependence of the Hessian matrix on the dynamical (field) variables and outputs the com¬ 
plete set of algebraically independent constraints. The Thomas decomposition splits the system 
of Euler-Lagrange equations into a set of involutive subsystems and automatically determines 
the complete set of constraints for each output subsystem which has a certain rank deficiency 
of the Hessian. Moreover, each output subsystem admits a formally well-posed Cauchy prob¬ 
lem 001 [i01 in terms of the output constraints. In contrast to other types of characteristic 
decompositions of (radical) differential ideals, a Thomas decomposition consists of differential 
systems whose solution sets are pairwise disjoint. This is accomplished by the use of inequa¬ 
tions (cf. Examr)le l4.6l) . allowing an irredundant description of the cases exhibiting different sets 
of constraints. 

This paper is organized as follows. Section |2] contains a short description of our research 
objects. In Section [2 we outline the standard approach I!]] to the computation of Lagrangian 
constraints. Our approach to this problem is based on the general algorithmic technique of dif¬ 
ferential Thomas decomposition jol llOl] outlined in Section |4] In Section |5] we adopt this tech¬ 
nique to the application to constrained dynamics and give a rigorous algebraic description of this 
adaptation. In Section |6] we illustrate its application by a number of examples including both 
field-theoretic and dynamical system models. Some concluding remarks are given in Section|7] 


2. Preliminaries 


In the framework of their field-theoretical description most of the fundamental laws of physics 
can be understood via the action 00 

where i] is a Lagrangian density depending on the field variables 

ip" - ... ,Xn,t) , a 6 {1,..., m), 
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and their first-order partial derivatives. In terms of its density X the Lagrangian L is expressed 
as 

The principle of least action (Hamilton’s variational principle) implies the (field) Euler-Lagrange 
equations 

p Ci p 

xq = t. ( 1 ) 

d{dx^<p‘‘) dtp^ 

Here and in what follows we use Einstein’s summation convention, i.e., summation over all 
possible values of the repeated indices occurring in a single term is assumed. 

Eor dynamical systems the action and the Euler-Lagrange equations read Qlli 


-I 


S= \ dtL(q“,q“), a e m] 


dt \ dq'^ j dq^ 


where q° are (generalized) coordinates and s dcf jdt are their velocities. 
A Lagrangian is regular if its Hessian matrix or Hessian defined as 


Hij = 


d^L 
dq\ dqi 
d^£ 
dip' d(fi 


(Dynamical System), 
(Eield-Theoretic Model), 


is invertible, and singular otherwise. 


( 2 ) 


(3) 


Remark 2.1. If a Lagrangian L possesses a continuous group of local gauge transformations, 
then its Hessian Q has a non-trivial nullspace. This fact follows from the generalized Bianchi 
identity (cf. ( 11 ]). Thus, gauge invariant theories are unavoidably singular. 


Let be a differential field of characteristic zero. By definition, K is a field containing the 
rational numbers Q with commuting derivations So, 6\, ..., 6„, i.e., additive maps Si : K ^ K 
satisfying the Leibniz rule 


Si{k\ k 2 ) - 6i{ki) k 2 + k\ Si(k 2 ) for all A:i, ^2 6 K, and S, o Sj - Sj o Si for all i, j. 

Moreover, let R„ - be the partial differential polynomial ring with differential 

indeterminates ff”, endowed with the set A = {5,, dx^,d,^^,..., dxj of commuting deriva¬ 

tions. In other words, is the polynomial ring 

K[d‘,di\ . .■d^x„p'‘ I Z>o,^ 6 {1,.. .,m)] 

with infinitely many indeterminates, where the action of dt, dx,, dx 2 , ■ ■ ■, dx„ on an indeterminate 
increases the corresponding exponent i, ji, ..., j„ by one and on K coincides with the action of 
So, Si, ..., 6 „. We identify ip' with ... d^yp' and also indicate differentiation by subscripts; 
e.g., pltx, ~ djdxtP^- We shall denote by Rq the ordinary differential polynomial ring with 
differential indeterminates ip^, ..., (p"' and A = [dt). 

In our study of the Euler-Lagrange equation systems ([TJ and (l2]i we shall restrict our consid¬ 
eration to field-theoretical and mechanical models such that 


£ e R„ and L e Rq, respectively. 
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It should be noted that this restriction is consistent with most of fundamental physical field theo¬ 
ries Qi] whose (singular) Lagrangian densities are differential polynomials in the field variables 
and their first-order partial derivatives over a field of constants. 

3. Linear algebra based approach 

In terms of the Hessian (O the set of Euler-Lagrange equations in ([T]l and (|2]l can be written 
as 

f Hi j qlj + Pi (Dynamical System), 

(5) 

Hi j (fil + Pi (Field-Theoretic Model), 

where the differential polynomials Hij and f*, contain derivations in t of order at most 1. 

If the Hessian H is singular, then there are differential consequences of equations in (|5]l that 
do not contain the second order derivatives in t. These consequences, called constraints, are of 
mincipal importance in the study of singular mechanical and field-theoretical models (cf. ifltl- 
iSil). Given a Lagrangian, the standard approach used in physics and mechanics to compute 
Lagrangian constraints (cf. Ill]) can be formulated as the following procedure: 

Computation of constraints based on linear algebra 


Step 1 Compute the Hessian H in accordance to (O. Then derive the set E of Euler- 
Lagrange equations (|5]l of cardinality m := lEI and let C := {). 

Step 2 Compute the rank r of the Hessian taking into account equations in E. 

Step 3 If r = m, then go to Step 6. Otherwise, go to the next step. 

Step 4 Compute a basis V of the nullspace of //, set up 

C:={P,y' la-1,...,|y|) 

and enlarge the equation set 

£:-£U{c = 0|c6C\{0)). 

Step 5 Set m r and go to Step 2. 

Step 6 Return C. 


The elements in C that have differential order at most 1 are called Lagrangian (cf., for exam¬ 
ple, fflEl), and the set of Lagrangian constraints is formed by a largest functionally indepen¬ 
dent subset. 

The above procedure suffers from the following algorithmic drawbacks. 

(a) The entries of the Hessian are ordinary or partial differential polynomials, i.e., we have 
Hi j 6 Rq or Hi j 6 Rn (cf. Section|2]). If Gaussian elimination is applied in Steps 2 and 4, 
the computation is performed in the field of fractions of Rq or in general. Division by 
a non-zero differential polynomial which vanishes on the solution set of the given system 
has to be prevented. Such differential polynomials are not detected by the above procedure 
(cf. also|(c)|below). 
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On the other hand, the rings Rq and R„ are not principal ideal domains. For this reason 
the Hessian H does not admit a Smith normal form in general, which would provide an 
algorithmic way to compute rank(//) in Step 2. In the case of a multivariate polynomial 
matrix its Smith normal form can be defined and computed only in very special cases 
(cf.lQ and the references therein) that are not related to the Hessian of the general form 
O). Similar remarks apply to the computation of a basis of the nullspace in Step 4 (cf., for 
example, El and its bibliography). 

(b) Generally, and we illustrate this fact by examples in Section|6] the rank of the Hessian may 
vary from one solution subspace of the Euler-Lagrange equations O or (|2]) to another one. 
Accordingly, the set of Lagrangian constraints may depend on a solution domain. 

(c) Computation of rank(i/) in Step 2 assumes simplification (reduction) of the Hessian mod¬ 
ulo (solution of) the system of Euler-Lagrange equations and its extension, respectively. 
Algorithmically, such a simplification can be performed by means of a characteristic de¬ 
composition 111 ill of the radical differential ideal generated by the left hand sides of the 
Euler-Lagrange equations. Thus, the procedure has to be extended with such a kind of 
decomposition. 

(d) The output set C of constraints has to be further processed to extract the set of Lagrangian 
constraints. 

Below we show that differential Thomas decomposition 10 E3] applied to the Euler-Lagrange 
equations Q provides a fully algorithmic way to compute an algebraically independent set of 
Lagrangian constraints that takes into account the dependence of the rank of the Hessian on the 
solution domain. 


4. Thomas decomposition and simple systems 


We sketch the idea of Thomas decomposition of systems of polynomial ordinary or partial 
differential equations (ODEs or PDEs) and its computation. This method traces back to work by 
J. M. Thomas in the 1930s 116]. Lor more details, we also refer to 0E3. 

Let 

pi = 0, ..., Ps -0, qi ..., qh^O {s,he Z>o) (6) 


be a system of polynomial ordinary or partial differential equations and inequations, i.e., the left 
hand sides p\, ..., ps, q\, ■ ■., qh elements of the differential polynomial ring Rq or /?„ over 
the differential field K of characteristic zero. 

A fundamental problem is to determine all power series solutions of ((h]). Writing the un¬ 
known functions ip^, ..., (f"' (of either f or of f, xi, ..., x„) as power series with indeterminate 
coefficients, we are led to substitute this ansatz into ® and compare coefficients. However, 
this naive procedure does not take into account conditions on the Taylor coefficients of i/j', ..., 
ip'" which arise as non-trivial consequences by differentiating the equations in (l6ll, taking linear 
combinations of their left hand sides, etc. Therefore, system (O needs to be transformed into an 
equivalent system which incorporates all integrability conditions and which is in this sense for¬ 
mally integrable (cf. 0). Lor nonlinear systems of PDEs certain case distinctions are necessary 
in general, resulting in a finite family of formally integrable PDE systems (e.g., a differential 
Thomas decomposition) such that the union of their solution sets equals the solution set of the 
original PDE system. Differential Thomas decomposition has various applications, allowing, 
e.g., to decide membership to the radical differential ideal which is generated by the given PDE 
system or to solve certain differential elimination problems M- 
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In order to identify the most significant variable in each non-constant (differential) polyno¬ 
mial, a total ordering > on the set of (partial) derivatives of all orders of all unknown functions 
ip'” is fixed beforehand. It is assumed to respect differentiation and that differentiation 
increases the significance of terms. The ordering > is referred to as ranking. 

We distinguish two stages of the computation of a differential Thomas decomposition, which 
may be intertwined in practice: its algebraic part and its differential part. The algebraic part of 
Thomas’ algorithm eliminates variables (representing unknown functions and their derivatives) 
in equations and inequations by applying polynomial division and ensures square-freeness of 
all left hand sides. For these reductions the corresponding left hand sides are considered as 
univariate polynomials in their most significant variables, which are commonly referred to as 
their leaders, with coefficients that are polynomials in lower ranked variables. 

More precisely, if two distinct equations p, = 0 and pj - 0 have the same leader v and 
deg^(p,) > degj,(py), then there exist (differential) polynomials ci, C 2 , and r such that 

c\Pi-C2Pj^r (7) 

and r is zero or its leader is ranked lower than v or its leader is v and deg^(r) < deg^(pj). This 
pseudo-division eliminates the highest power of v in p,, where ci may be chosen as a suitable 
power of the initial init(pj) of pj, i.e., the coefficient of the highest power of v in pj. Replacing 
the equation p, = 0 with the equation r = 0 in the system requires that c\ does not vanish 
on the solution set of the system, if the solution set of the system is to be maintained by this 
operation. For this reason, Thomas’ algorithm possibly splits the system into two systems which 
are subsequently treated in the same way as the original system, which is discarded. These 
two systems are defined by inserting the inequation init(pj) + 0 and the equation init(pj) = 0, 
respectively. 

For each two distinct inequations qi + 0 and qj + 0 with the same leader the least common 
multiple (1cm) r of qi and qj is computed and the two inequations are replaced with r 0. 
This process involves pseudo-divisions as above and also distinguishes cases of vanishing or 
non-vanishing initials. 

For each pair p, = 0, qj + 0 of an equation and an inequation with the same leader the 
greatest common divisor r of p, and qj is computed, where polynomials are considered again as 
univariate polynomials in their leaders and it is required that initials of pseudo-divisors do not 
vanish on the solution set of the system. If p, divides qj, then the system is inconsistent and is 
discarded. If r e /f \ {0), then the inequation qj 0 is removed from the system. Otherwise, the 
equation p, = 0 is replaced with the equation p,7r = 0. 

Finally, further reductions and splittings of systems may be necessary to ensure square- 
freeness of all left hand sides (as univariate polynomials), i.e., non-vanishing of their discrimi¬ 
nants. Since K is of characteristic zero, square-freeness of a univariate polynomial is character¬ 
ized by the condition that the polynomial and its derivative have no root in common. Adhering 
to the recursive representation of polynomials, a computation of the greatest common divisor 
of a left hand side and its derivative with case distinctions as above allows to determine the 
square-free part. 

The above procedure is a variant of Euclid’s algorithm for multivariate polynomials, termi¬ 
nating after finitely many steps. If the fact that the variables occurring in p, and qj represent 
unknown functions and their derivatives is neglected, then the result is a finite collection of alge¬ 
braic systems with the following property. 
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Definition 4.1. An algebraic system S as in (|6]l is said to be simple (with respect to >), if the 
following three conditions are satisfied, where is the subsystem of S consisting of those 
equations and inequations whose leader is ranked lower than the variable v. 

(a) All Pi and all qj are non-constant polynomials. 

(b) The leaders of all p, = 0 and qj + 0 are pairwise distinct. 

(c) If V is the leader of p, = 0 or qj + 0, then neither the initial nor the discriminant of that 
equation or inequation has a solution (over the complex numbers) in common with the 
subsystem 5 

Remark 4.2. A simple algebraic system S can in principle be solved recursively over the com¬ 
plex numbers by extending the coordinates of a solution of S with a root of the equation with 
leader v, if present, or with a complex number satisfying the inequation with leader v, or with an 
arbitrary complex number if v is not the leader of any equation or inequation in S. 

Moreover, conditions |(a)| and |(b)| imply that the left hand sides pi, ps, qi, qi, in a 
simple algebraic system are algebraically independent over K, i.e., the zero polynomial is the 
only polynomial with s+h indeterminates and coefficients in K which vanishes under substitution 
ofpi, ...,p„qi, ...,qh. 

The differential part of Thomas’ algorithm works towards formal integrability. For each 
two distinct equations p; = 0 and pj - 0 whose leaders d'‘(p‘‘ and d‘(p‘‘ 

gf are derivatives of the same unknown function the algorithm considers the 
derivatives of p, and pj which both have leader lcm(5*, d^') (ff. In an appropriate linear combina¬ 
tion of these two derivatives the highest power of this leader is eliminated. Such a linear com¬ 
bination is referred to as a cross-derivative of p, and pj. The algorithm checks whether pseudo¬ 
division of this cross-derivative modulo (derivatives of) the equations of the system yields the 
zero polynomial. Non-zero remainders resulting from this reduction are added as new equations 
to the system, and the augmented system is treated by the algebraic part of Thomas’ algorithm 
again. A method which traces back to M. Janet lO allows an efficient organization of all cross¬ 
derivatives. By distinguishing admissible and non-admissible derivations for each equation, each 
cross-derivative is considered only once. 

The differential part performs pseudo-divisions as in O, which however involve derivatives 
of equations of the system in general. If the equation p, = 0 is replaced with the equation r = 0, 
then we require that ci does not vanish on the solution set of the system. If a proper derivative 
of Pj is subtracted, ci may be chosen as a suitable power of the initial of that derivative, which 
is the separant of pj, namely, the partial derivative of pj with respect to its leader. Since the 
discriminant of a differential polynomial (with respect to its leader) is essentially the resultant of 
the polynomial and its separant, the algebraic part of Thomas’ algorithm ensures that separants 
do not vanish on the solution set of the system. 

The result of Thomas’ algorithm is a finite collection of differential systems with the follow¬ 
ing property. 

Definition 4.3. A differential system S as in (|6]l is said to be simple (with respect to >), if the 
following three conditions are satisfied. 

(a) The system S is simple as an algebraic system (with respect to >, cf. Def. l4.1l) . 

(b) The cross-derivative of each pair of distinct equations whose leaders involve the same 
unknown function reduces to zero modulo the equations of the system and their derivatives. 
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(c) The left hand side of every inequation is reduced modulo the left hand sides of the 
equations, in the sense that no pseudo-division of qj modulo any p, is possible. 

Remark 4.4. Note that each splitting of a system S into two systems S \ and S 2 performed by 
Thomas’ algorithm is arranged so as to ensure that the solution set of S is the disjoint union of 
the solution sets of S \ and S 2 - 

Definition 4.5. Let 5 be a differential system as in A hnite collection S i, ..., Sk of sim¬ 
ple differential systems whose solution sets form a partition of the solution set of S is called a 
differential Thomas decomposition of S (with respect to >). 

A Thomas decomposition is not uniquely determined and its computation depends on the 
choice of the ranking > and many other choices. Among the possible improvements of the above 
outline of Thomas’ algorithm is the option to factorize the left hand sides of equations and to 
split a system according to a non-trivial factorization. 

The algorithmic details of the Thomas decomposition method have been worked out in 
(cf. also ilOr i. Implementations have been developed by T. Bachler and M. Lange-Hegermann as 
Maple packages AlgebraicThomas and Dif f erentialThomas and are freely available M- 

We illustrate the Thomas decomposition technique on a simple example. 

Example 4.6. Let us consider the PDE given by the following Hessian determinant; 



Ux,x 

^x,y 

^X,z 

det 

^x,y 

^y,y 

^y,z 


s ^x,z 

^y,z 

^z,z > 


( 8 ) 


The Thomas decomposition method is applied over the differential polynomial ring K{u} with 
one differential indeterminate u and commuting derivations dx, dy, d^, where AT is a differential 
field of characteristic zero, e.g., K - Qso that the restrictions of dx, dy, d^ to K are all zero. We 
choose the degree-reverse lexicographical ranking satisfying dxU > dyU > d^u, i.e.. 


dl'dP/df u > 5*'5^^5*^ u 


Q.\ + 0-2 "I" ^3 ^ + ^2 "I" ^3 

(Qi Q 2 + ^^3 — b\ b 2 + b^ and 
{a\,a 2 ^ai) {b\^b 2 -,bi) and 

ai < bi for i - max {j \ aj ^ bj }). 


Then the initial of the PDE (O is u^,z‘^y,y ~ “yz- We distinguish the cases whether this initial 
vanishes or not. In the first case we obtain the following simple differential system under the 
assumption that the initial ^ of the initial of (O does not vanish (no integrability condition 
needs to be checked for the single equation); 


(Ti) 


^y,y ^y ^x,x ^ ^y^z ^x,z ^x,y ^z,z v ^y^y 


- 

0 

^z,z ^y,y ~ 


+ 

0 


^z,z 

+ 

0 


Leaders of differential polynomials are underlined where not obvious. If j. vanishes in the first 
case, a simplihcation yields the second simple differential system 


(Ti) 


~^y,Z "t" 2 Uy z tix^l tlx,y tly^y ^ ~ 0 , 

tly z ^ 0 , 

M77 = 0 
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after pseudo-reduction of a single cross-derivative to zero modulo the system. If the initial of (l8]l 
vanishes, we again distinguish the cases whether vanishes or not. If not, the pseudo-division 

(o — 2 \ , 1 ( — ^ \ — I _ 

^z,z ^y,z ^x,z ^y,z ^z,z ^x,y ^y,y ^x,z) ^x,z \^z,z ^y,y ^y,z) ~ y^z,z ^x,y ^y,z^x,z) 

allows to replace the PDE ([8]l with the expression in the bracket on the right hand side (ensuring 
square-freeness). We obtain the simple differential system 


(Ti) 


^z,z ^x,y ^y,z ^x,z 

.,2 


0 , 

0 , 


^Z,Z ^ 0 


after confirming that the integrability condition for the two equations is satisfied by applying 
pseudo-reduction modulo the system. Finally, if the initial of the PDE ([8]l vanishes, then = 0 
implies Mj, ^ = 0, and (O reduces to Uy^ — Distinguishing the cases whether Uy^y vanishes 
or not, we obtain the simple differential systems 


and 


(7’4) 


{Ts) 


M.l.Z = 0 , 

^y,z ~ 0 5 

^z,z ~ 0 5 

y ^ 0 

My y - 0 , 

* Ity^y^ - 0 , 

^z,z ~ 0 . 


Thus, we have the following tree of case distinctions (suppressing branches leading to inconsis¬ 
tent systems): 



As opposed to the Thomas decomposition (Ti)-(Ts) with pairwise disjoint solution sets, e.g., 
the package Dif f erentialAlgebrain Maple 17 computes, using the Rosenfeld-Grobner algo¬ 
rithm, the following characteristic decomposition for ® : 

i iyz,z ^y,y 2 Uy y^ ^x,z ^x,y ^z,z ^x,y ^yy ^x,z ~ ^ ’ 

I j My y — My^^ 0 
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'Z,Z 


0 , 

0 


(=(7’3)) 




2 

^y,y ~ ^y,z 
^z,z 


+ 


U-X,Z — 0 > 

* lly^^ — 0 , 

^Z,Z ~ 0 


^y,y — 

< Uy,, = 0, (=(r5)) 

^Z,Z ~ 0 ■ 


The solution sets of these differential systems are not disjoint. The third and fourth systems have 
solutions /i(x) + / 2 (x)y + cz in common, where /i and /2 are arbitrary analytic functions and c 
is an arbitrary constant. 


Remark 4.7. Due to its disjointness the Thomas decomposition partitions the solution space 
of the Euler-Lagrange equations ([T]i or (|2]). In doing so, it may occur that the field-theoretic or 
mechanical model, which is globally non-singular, is singular in some domains of the solution 
space. In those domains there are hidden Lagrangian constraints that are revealed by the Thomas 
decomposition (cf. Section |6]for examples). In terms of these local constraints one can pose an 
initial value (Cauchy) problem providing existence and uniqueness of local analytic solutions 
T2]. The local analyticity provides a smooth dependence of the solution on the initial data 
and thereby makes the initial-value problem well-posed (cf. lll9n i. It should be noted that the 
well-posedness of the initial value problem has been proven for an involutive simple differential 
system ||3l related to the domain under consideration. For another type of characteristic decom¬ 
position which is not disjoint a proof of the well-posedness for Cauchy problem is not known in 
the literature. 


Remark 4.8. In some applications one might be interested only in one branch of the tree of case 
distinctions. For instance, if the differential ideal which is generated by the left hand sides of the 
given system of PDFs is prime, the resultin g Th omas decomposition contains a simple system 
which is most generic in a precise sense (cf. lIlOl Subsection 2.2.3]). Clearly, if only the generic 
simple system is required, the computation can be restricted to that branch. 


5. Adaptation of Thomas decomposition to constrained dynamics 

Hereafter, we shall write the underlying ring of differential polynomials in (HJi without sub¬ 
script, since it will be clear from the context what ring is meant, partial or ordinary. We fix a 
ranking > onR which satisfies 

dt ip'^ > • • ■ 5*" I/)* for all ki,..., k„ e Z>o, a,b & [I,... ,m]. (9) 

The leader of any non-constant differential polynomial p in R, i.e., the differential indeterminate 
or derivative 5/(9^', ■ ■ - ipf which is maximal with respect to > among those occurring in p, is 
denoted by ld(p). 
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Remark 5.1. Condition (|3) implies that for all i, j, ki,... 6 Z>o and for all a,b e 

{1,..., m) we have 

j>k ^ 

The claim follows from applying to the left hand side and the right hand side of 

Definition 5.2. Let ip e ,..., be minimal with respect to >. A non-constant differential 
polynomial p e Ris said to be a -condition if 

id(p) < d^iif 


holds. 

Definition 5.3. Given Euler-Lagrange equations ([TJ or (l2]l, their differential consequence p - 0 
is called generalized Lagrangian constraint if p is a 5^i/'-condition. If a constraint holds for all 
solutions of O or © it is called global, and local otherwise. A non-singular Lagrangian system 
with local constraints is locally singular. 

Remark 5.4. For singular models of mechanics, generalized Lagrangian constraints are just 
Lagrangian ones. For singular (or locally singular) field-theoretical models, Lagrangian con¬ 
straints, as we underlined in Section [3] being of the first differential order, form a subset of the 
set of generalized Lagrangian (or local Lagrangian) constraints. 

We consider the subring 

^ordw <2 := K[d‘,di[ ■ \ i 6 {0, 6 Z>o,ke m]] 

= K[<f^, 5,,/,..., 5,/, ...\ke{l,...,m}] 

of as a differential ring with derivations ■■ ■, dx„ and the ranking induced by >. 

Remark 5.5. Every S^i/'-condition is an element of f?ord( 0 < 2 - 

Remark 5.6. Let / be a differential ideal of R, i.e., an ideal of R which is closed under all deriva¬ 
tions of R. The set of 5^i/r-conditions in I is closed under taking /?ord(f)< 2 -linear combinations and 
applying an arbitrary number of the derivations dxi, ■ ■ ■, dx„ (cf. Remark lSTTI) . In other words, the 
set of 5^t^-conditions in / is a differential ideal of Rord(t)< 2 - 

We denote by 5" the set of left hand sides of equations in a differential system S. The 
saturation of a differential ideal £ of R by a differential polynomial q e R is defined as the 
differential ideal E : q°° :- {p e R \ p e E for some k e Z>o ). 

Proposition 5.7. Let S be a simple differential system over R with respect to a ranking satisfying 
m and let f/ 6 [ip^, ..., ip”') be minimal with respect to >. Let E be the differential ideal of R 
which is generated by S'" and let q be the product of the initials and separants of all elements of 
S". Moreover, let 

S ord(f)<2 := {p 6 5 = I ld(p) < djif }, 

Eoid(t )<2 the differential ideal o//?oid(f )<2 which is generated by S ord( 0 < 2 . ond qoid{t )<2 the product of 
the initials and separants of all elements ofSoid{t)< 2 - Then the set of djf-conditions iniE : q°° 
is the radical differential ideal £’oid( 0<2 : (^ord( 0 < 2 )“ ofRoid(t)< 2 - 
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Proof. Let p be a ^^^/'-condition in I. Since 5 is a simple differential system, repeated differential 
pseudo-reduction of p modulo S ~ yields zero (cf. lH^ Prop. 2.2.50]). Since we have ld(/7) < 
only elements of 5ord(f)<2 are chosen as pseudo-divisors. Their initials and separants are elements 
of Rori{t )<2 and the coefficients which are used for the pseudo-reductions are in /?oid( 0<2 as well. 
This shows that the set of cl^^/r-conditions in I is contained in £’ord( 0<2 : (^oid( 0 < 2 )“- Conversely, 
every element of £’ord(f )<2 : (^ord(f)< 2 )°° is clearly a cl^i/r-condition in I. □ 

We recall that the radical of a differential ideal I of R 'k defined as the differential ideal 
V7 := { p e I 6 / for some k eN}. 


Theorem 5.8. Let 


pi=0, ..., ps=0, qi^O, ..., qh^O (10) 

be a (not necessarily simple) system of polynomial partial differential equations and inequations, 
where p\, p,, qi, , qh & R, s, h & Z>o. We denote by E the differential ideal ofR which is 

generated by pi, ..., Ps and by T the set of if-conditions in sjE : q°°, where q q\ ■ ... ■ qf,. 
Let 

be a differential Thomas decomposition ofm with respect to a ranking satisfying m and let 
if 6 {(^*,..., ip"'} be minimal with respect to >. For i — , r, let 

j-ii) .— p(0 . (W) 

• ^ord(f)<2 • ''Vord(0<2^ ’ 

where differential ideal o//?ord(o <2 which is generated by 

and product of the initials and separants of all elements o/‘5ord(f)<2' have 

T = r*') n... n 

Proof. We have 

yjE : q^ = (£^1) : (^<^’)“) n ... n (E^’''> : (q^''^D, 

where E^'^ is the differential ideal of R which is generated by (S^‘))= and is the product of the 
initials and separants of all elements of where i e {1,..., r) (cf. llOl Prop. 2.2.72]). From 

Pror)osition l5.7l we conclude that, for each i e {1,..., r], the set of cl^i/^-conditions in ; (q^‘^)^ 
is equal to T^‘\ □ 

In the examples of the next section it is a very simple task to trace the origin of the djif- 
conditions arising in the differential Thomas decomposition. 
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6. Examples 


In this section we apply the Thomas decomposition to four models. Two of these models are 
field-theoretic and the other two are dynamical systems. 


Sample 6.1. The first model is the (1H- l)-dimensional version of the chiral Schwinger model 

In the framework of Hamiltonian formalism the last model was studied, for example, 




221} . We consider the Lagrangian density in the form ( Il22l] . Equation 19) 




i (d,Ao - -I- i (d,^f - ^ + e id,(p)Ao + e 4>{d,Ai) 

+ e (5,0) (Ao-Ai)+^ae^ (Ag - Aj). 


Here, e, a are parameters, f, x are the independent variables and i/j' = Aq, - Ai, (p? - (p are 
the dependent variables. 

Let > be the ranking satisfying (|3l and such that 


^ < A\ < Ao < dx<p < dxA\ < 5jAo < d\tp < ... 

< 5,0 < 5,Ai < 5,Ao < 5,5,v0 -< dtdxA\ < dtdxAo < dtd\(p < ■ ■ ■ 

The Euler-Lagrange equations read 


5^Ao - 5,5xA 1 - e (dttp -i- 5x0) - a e^Ao 

= 0, 


5,5xAo - e (dt4> + 5x0) - 5xAi - a Ai 

= 0, 

( 11 ) 

5^0 -H e (5,Ao - 5,Ai) - 5^0 -H e (5xAo - 5xAi) 

= 0. 



The Hessian (l3]l of this field-theoretic model is H - diag(l,0,1) and, hence, the model is 
singular. The Euler-Lagrange equations (fTTT) are linear. In this case Thomas’ algorithm applied 
to the system performs its completion to involution and outputs the Janet involutive form of (fTTT i 
(cf. l^liOl]). Denote the first and second equation in (fTTT i by P and Q, respectively. Then the first 
step of the completion procedure consists in detection of the (underlined) leaders, calculation of 
the cross-derivative 

dxP-d,Q^0, 

and then elimination (reduction) of terms in the left hand side modulo equations in (fTTT i. as 
explained in Section|4] As a result of the completion, Thomas’ algorithm yields 

(1 — a) 5,Ao -I- (1 -I- fl) dxAfs — 5,Ai — dxAi = 0, 

(1 -I- a) (5?Ai - 5^Ai) - e (2-H a) (5,0 -I- 5i0) - a (Ao -I- Ai) - Ai = 0, 

i , . (12) 

(a -I- l)(5,5^Ai - 5^Ao) - e (5,0 -I- 5x0) - a Ai = 0, 

5^0 - 5^0 - ea (5.,Ao - 5,Ai) = 0. 

The first equation in (fT2t (marked in boldface) is a Lagrangian constraint. 
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Example 6.2. To show a field-theoretic situation when the differential Thomas decomposition 
computes local (generalized) Lagrangian constraints (cf. Definition 15.31) for a (globally) non¬ 
singular (but locally singular) model, consider the following Lagrangian density: 




-U, + UV,Wx + W, (u, + Vx) ■ 


(13) 


The independent variables are f, x and the dependent variables are - u, - v, (p^ - w. We 
choose the ranking > satisfying (|3 such that 


W < V < U < Wx < Vx < Ux < Wx,x < ■ ■ ■ < Wt < Vt < u, < Wt,x < Vt,x < Ut^x < Wt,x,x < ■ ■ ■ 


The density (fTlT l generates the Euler-Lagrange equations 


Ut,t + Wt,! -V^Wx - 0 , 

2uWx Vt^t + 2utVtWx + 2uvt Wi,x + Wt^x — 0 , 

Ut,t + Vt,x + 2uVt Vt^x + Uxvj - 0 . 


(14) 


Let A, B, C denote the first, second and third equation, respectively. Considering the (underlined) 
leaders of these equations, Thomas’ algorithm first replaces C with C - A. This corresponds to a 
row reduction of the Hessian: 




1 0 1 ’ 
0 2 uwx 0 

1 0 0 , 


The leaders of the Euler-Lagrange equations are then Vt,t, Wt,t. Since these differential poly¬ 
nomials involve pairwise distinct dependent variables, formal integrability of the system requires 
no further check (i.e., technically speaking, all derivations are admissible in the sense of Janet 
division). Only the possible vanishing of the initial of the second equation has be taken into 
account (compare also with the entry at position (2,2) of the Hessian i/*''). 

The field model (fT4l l is not singular. Eor this reason the generic case leads to a simple system 
not containing any constraints (cf. (Ti) below). We conclude, using Theorem l5.8l that among the 
consequences of the given PDE system (fTTt there exist no (generalized) Lagrangian constraints. 


There are precisely two possibilities for the initial of the second equation to vanish: Wx = 0 
or u - 0. All local generalized Lagrangian constraints for these more special consequences of 
the given PDE system (fT4l) arise from the imposed equations or their integrability conditions ||3l. 

In the first case (i.e., Wx - 0), the integrability condition for the non-zero equations in (fT4l) 


Wx - 0 , + Uxvj + 2uVt V,^x - Wt,t — 0 


is 

(2 M V, -I- 1) Vt,x,x + 2 M + 4 Ux Vt V,,x + Mx,x vj - 0 . 

We obtain three simple systems distinguishing the cases whether the initial 2uvt + 1 of the 
last equation vanishes (cf. (T 4 ) below), or has vanishing initial u (cf. (Tj,) below), or has non¬ 
vanishing initial u (cf. {T 2 ) below). 

In the second case (i.e., u-^,yvx + 0), the integrability condition for 

Wx = 0 , 

~ 14 


w,,x = 0 






is 

2 vtWxV^ = 0. 

We obtain two simple systems distinguishing the cases whether v, vanishes or not (cf. (Ts) and 
(Te) below, respectively). 

Summarizing we obtain the following tree of case distinctions 



where (Ti), (T(,) are the simple differential systems of the following Thomas decomposi¬ 

tion of (fT4l l and where for every simple subsystem we mark its local generalized Lagrangian 
constraints in boldface: 


(Ti) 


U,j + (2uVt+l)v,^:, + UjcV^ 
2uWx Vt^t + (2uvt + 1 ) Wt,x + 2.WxVt u, 
< (2uVt + 1) V/,x - Wij + (Ux + Wx) vj 

Wx 

u 


= 0 , 

= 0 , 

= 0 , 

* 0 , 

0 


{T 2 ) 


U,^, + {2 UVt + 1) Vf_x + Ux vj 

= 

0, 

U + 4 UxVt Vt,x + vj Ux,x 

= 

0, 

(2uVt+ 1) V,^x + Ux v] - Wtf 

= 

0, 

Wx 

= 

0, 

u 


0, 

2uvt + 1 


0 


(Ti) 


(Ta) 


M = 0, 

Vt^ = 0, 

- 

Wi^t - 0 , 

Wx = Q 

Utj - 0, 

Ux = 0, 

2h -I-1 = 0, 

Wt,t - 0 , 

Wx = 0, 

u + 0 
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(Ts) 


(Te) 


< 


H = 0, 

Vf = 0, 

< w,_, = 0, 

= 0, 

Wx ^ 0 

u = 0, 

V<,* = 0, 

Wt,t -WxV^ = 0, 

= 0, 

Wx,x = 0 , 

Wx + 0 , 

V, + 0 . 


Since the differential Thomas decomposition is disjoint (cf. Section|4li, the simple differential 
subsystems (Ti),.. (Te) (cf. Definition |43]) partition the solution space of the Euler-Lagrange 
equations (0, and reveal local Lagrangian constraints. 


Example 6.3. As an example of singular dynamical system we consider Lagrangian taken from 
( |@1, Equation 8.1) 


^ - ^2 (^l)? + + 2 ^1 ^2 {qi)t {q2)t + + ?2 • 

The independent variable is t and the dependent variables are - q\,y^ - qi- We choose the 
ranking > such that 

q 2 <qi< {q 2 )i < {qi)t < (q 2 )t,t < (qi)t,t < ■■■ 

Then the Euler-Lagrange equations are given by 

( 4q2(q2)t(qi)t+ 2ql(qi\,+ 2 qiq 2 (q 2 )t,t-2qi = 0, 

1 2 - ( 15 ) 

[ 4qi(q2)t(qi)t + 2qf(q2)t,t + 2qiq2(q^-2q2 = 0 . 

Let A and B denote the first and second equation, respectively. According to the (underlined) 
leaders of these equations, Thomas’ algorithm first performs a pseudo-reduction of B modulo A; 
more precisely, B is replaced with the remainder 

q2 B-qxA^2 {q\-ql)^0. (16) 

This computation corresponds to a row reduction of the Hessian: 

^(1)^/ '^^2 2^i^ 2 'l 

\ 2 ^i ^2 2 .q\ I ' 

We obtain a zero row (indeed, the determinant of vanishes), which is reflected by the fact 
that all terms involving differentiation order 2 in (fThl l cancel. Hence, the model (fTSl l is singular. 

The solution set of the system does not change when B is replaced with the remainder in 
(fThl l if the coefficient q 2 in (fThl l does not vanish on the solution set. Eor this reason, Thomas’ 
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algorithm splits the original system (Eli into one system incorporating the condition q 2 + 0, 
where B is replaced with q\-q\ - 0, and a complementary system containing the new condition 
q 2 - 0. By taking the factorization q^ - q^ - iqi - q 2 ){q\ + qi) into account, the first system is 
split again into two complementary systems iT\) and {T 2 ). For both systems, Thomas’ algorithm 
reduces A modulo q\ - q 2 - 0 o^: q\ + q 2 — 0, respectively. Both remainders are divided by the 
non-vanishing factor 2 q 2 . Note that for ODE systems no formal integrability check is necessary, 
and the initials of the resulting left hand sides in (Ti) and {T 2 ) do not vanish. By reducing A and 
B modulo q '2 = 0 in the remaining case, we obtain the third simple differential system {T^) of the 
Thomas decomposition. The Lagrangian constraints we find are qi - q 2 - 0 and qi + q 2 - 0 in 
(Ti) and {T 2 ), respectively (in boldface below). 

Thus, we have the following tree of case distinctions 


9i - 92 = 0 


(Fi) 


(El 



where (Ti), (T2), (T^) are the simple differential systems of the following Thomas decomposition 

ofdnii: 


(Ti) 

(T 2 ) 


2 q 2 (q 2 )t,t + 2 (q 2 )j - 1 = 0 , 

qi-qi = 0, 

qi * 0 

' 2 q2 {q2)i,i + 2 (q' 2 )? - 1 = 0 , 

^ qi + qi = 0, 

qi * Q 


(Ti) 



0 , 

0 . 


Note that the local Lagrangian constraints in the simple systems (T 1 ) and (T 2 ) can be combined 
in a single global constraint q^ — q^ = 0. 


Example 6.4. The double sombrero model ( iEbIi . Equation 7). Its Lagrangian is given by 

^ (^2)? + (qi)j - kf + ^qql - ^Aq'l. 

The independent variable is f, the dependent variables are = qi, = qi, and k, A, and fi are 
non-zero parameters. We choose the ranking > such that 


q2<q\< iq2)t < (qi)t < iq2)t,t < (qi)t,t < ■■■ 
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The Euler-Lagrange equations read 


2 ((gi)f {qi)t,t + q\ (q\)t{q2ft + q\ (q2)t(q2)t.t) (qi)t + 

{(qi)f + qi(q2)j-kj{(q^-qi(q2)j)-IJ.qi+Aql = 0 , 

2 ((^1)/ (qi)t,i + qiiqi)t (q2)j + q\ (?2)r {q2)t,) q\ {q2)t + 

{{qif;+ q\{q2)]-k){2qi{qx),{q2)t + q\{q2)t.) = 0 . 

The Hessian is 

( 1 ) q\(q 2 )^ + '^iq\)li -k ' 2 .q\(,q\)t{q 2 )t 

2 q\{qi),{q 2 ), 2 q\{q 2 )^ + {q\(,q 2 )^ + {qi)^ - k)q\ ^ 

and its determinant is 

det//^^> = q\ (q'j (^2)? + {q\f, - (3 q\ (^2)? + 3 {qif, - k'j . 

There are solutions of the Euler-Lagrange equations (fTTl i such that the Hessian does not vanish. 
Therefore, the model is not singular, and there are no global Lagrangian constraints. 

Let A and B denote the first and second equation in (fTTl l. respectively. We have 

init(A) = q\ (^ 2 )? + 3 {qif, - k , init(B) ^ 2 q] (qi), (q2), ■ 

Under the assumption that init(A) and its separant 6 {qi)t do not vanish on the solution set of the 
system, Thomas’ algorithm replaces B with the equation init(A) B - init(B) A = 0, whose leader 
is {q2)r,r and whose initial is 

qi {q\ (^ 2 )? + (?i)? - k) (3 q\ (^ 2 )? + 3 {qx)^ - kj . 

The coefficient of {q2)t,t in A is 2 q^ {qi)t (q2)t- Therefore, the algorithm subsequently replaces A 
with 

{q\ (q2)j + (qi)j - (3 ql (^ 2 )? + 3 (qi)^ - k^A + 2 qx (qx), (q2)t B - 0 , 

whose left hand side is (exactly) divisible by the initial q\iq 2 )j + 2iiqx)j - k of the previous 
equation A. Under the above assumption, we may divide by this expression and replace the 
left hand side with the quotient. Since the leaders iqx)t,t and (q 2 )t,t of the new equations A and 
B are derivatives of distinct unknown functions, no integrability check is necessary. Under the 
assumption the initials of these equations and the discriminant of these initials (which is a certain 
multiple of (q^ (^ 2 )? - k){3 q\ (^ 2 )? - k)) do not vanish, we obtain the simple differential system 
ijx) below. The different cases of vanishing initials or discriminants have to be treated and this 
yields four other simple differential systems {T 2 )-{Ts). Note also that two pairs of factors of the 
inequation in (Ti) become equal if {q2)t is specialized to zero, which is the reason for a case 
distinction leading to (7’2). We obtain the following tree of case distinctions 
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where we have only displayed cases leading to consistent systems and where (Ti), (T^) 

are the simple differential systems of the following Thomas decomposition of (fTTl) with local 
Lagrangian constraints marked in boldface; 


(Ti) 


- (q 2 )j + (qi)j - k) (3 q\ (^ 2 )? + 3(q'i)2 - k) (qi\, + 

3^1 iqi)^ (qi)t + qi ( 6^1 (qi)^ -^k{q 2 )j - Aq\+^j) (qi)^ + 

3 q\ {qift -^kq\ (q 2 )‘f - q\{i Aq\ - 3 Iiq\ - k^) (^ 2 )? + 

kAq\ — k^q\ - 0, 

[q\ iqi)^ + (?i)? - k) (3 q\ (^ 2 )? + 3 (g'l)? - k) qi (q 2 \t + 
6 (? 2 )r(g'l)f +4(3 q]{q 2 )j - 2 k)(^2)?(?l)? + 
{6q‘l(q2)^-^kq'^^(q2)'^-2{^Aq‘l-i2q]-k^'j(q2)l)(qi)i = 0, 
qi (q 2 )t {q\ (qift - k) {q\ (q 2 )j + (qi)j - k) 

{q\{q2)^+ 3{qx)^-l^(3q\{q2f, -l^[3q\(q2)^+ 3{qx)^-1^ + 0 . 


Note that the initials of the left hand sides of the above two equations divide det//4) 


(7’2) 


[3(q{)]-k'^{qi),j + Aq\-^iqi = 0, 

{qi), = 0 , 

q\{{q\)^-k)(3(qi)^-1^ + 0 


(T2) 


= 0 


(7’4) 


(7’5) 


^Aq\-^yiq\ + lk^ 

3 q\ iqi)] - k 

Aq\-n = 
kt iqz)] -kA = 


0 , 

0 


19 


0 , 

0 . 










Remark 6.5. We note that in our examples Thomas’ algorithm is applied to a system of non¬ 
linear differential equations which are linear in their leaders, i.e., in their highest derivatives 
with respect to the chosen ranking >. Therefore, one of the first steps in a computation of a 
differential Thomas decomposition for the given PDE system corresponds to the row reduction 
of the corresponding Hessian, which is part of the approach based on linear algebra recalled in 
Section |3] 


7. Conclusion 


In this paper we have shown that differential Thomas decomposition applied to the Euler- 
Lagrange equations of field-theoretical ([TJ or mechanical 0 m odels with polynomial Lagrangi- 
ans is a proper algorithmic tool, implemented in Maple llSt], for computing algebraically in¬ 
dependent Lagrangian constraints or generalized Lagrangian constraints, respectively. In doing 
so, the decomposition being disjoint generates a partition of the solution space of the Euler- 
Lagrange equations. In the case of a singular or locally singular model the complete set of the 
corresponding local and algebraically independent constraints is computed for every element in 
the partition. 

One can also apply another characteristic decomposition of the radical differential ideal gen¬ 
erated by the Euler-Lagrange equations, for example, that based on the Rosenfeld-Grobner algo¬ 
rithm and implemented in the dif f alg lib rar^ of Maple (cf. m and its bibliography). Maple 
has one more built-in splitting algorithm 12411 for nonlinear systems of differential equations 
(command rif simp, a part of the package DEtools), that can also be used for the detection of 
Lagrangian constraints. However, diffalg and rif simp do not yield disjoint decompositions 
in general, and the Lagrangian constraints that are inherent to different output subsystems may 
interfere. 

The ranking Q is a Riquier ranking ^ (cf. also ||3l)- Thus, by the Riquier existence 
theorems i27ll28|] . a simple differential subsystem with nonempty equation set provides the ex¬ 
istence and uniqueness of formal power series solution satisfying certain initial (Cauchy) data 
(cf. fiQ). If the simple system under consideration is singular, and hence contains (general¬ 
ized) Lagrangian constraints, then their presence is to be taken into account by the initial data, 
and the second-order derivatives in time (‘accelerations’) are uniquely defined. 

Differential Thomas decomposition can also be applied to detect and compute hidden con¬ 
straints in differential-algebraic (and also in partial differential-algebraic) equations (DAEs). The 
presence of hidden algebraic constraints, i.e., those algebraic constraints that are not explicitly 
given in the system, is the main obstacle in numerical solving of DAEs (cf., for example, ll2^ 1. 
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